Rasters.jl

RastersModule

Rasters

CI Codecov Aqua.jl Quality Assurance

Rasters.jl defines common types and methods for reading, writing and manipulating rasterized spatial data.

These currently include raster arrays like GeoTIFF and NetCDF, R grd files, multi-layered stacks, and multi-file series of arrays and stacks.

EarthEnv HabitatHeterogeneity layers trimmed to Australia

A RasterStack of EarthEnv HabitatHeterogeneity layers, trimmed to Australia and plotted with Plots.jl

Data-source abstraction

Rasters provides a standardised interface that allows many source data types to be used with identical syntax.

  • Scripts and packages building on Rasters.jl can treat Raster, RasterStack, and RasterSeries as black boxes.
    • The data could hold GeoTiff or NetCDF files, Arrays in memory or CuArrays on the GPU - they will all behave in the same way.
    • RasterStack can be backed by a Netcdf or HDF5 file, or a NamedTuple of Raster holding .tif files, or all Raster in memory.
    • Users do not have to deal with the specifics of spatial file types.
  • Projected lookups with Cylindrical projections can by indexed using other Cylindrical projections by setting the mappedcrs keyword on construction. You don't need to know the underlying projection, the conversion is handled automatically. This means lat/lon EPSG(4326) can be used seamlessly if you need that.

Named dimensions and index lookups

Rasters.jl extends DimensionalData.jl so that spatial data can be indexed using named dimensions like X, Y and Ti (time) and e.g. spatial coordinates.

Regions and points can be selected with a..b, At Near and Contains.

Dimensions can also be used in most Base and Statistics methods like mean and reduce where dims arguments are required. Much of the behaviour is covered in the DimensionalData docs.

See the docs for more details and examples for Rasters.jl.

source

Common Applications

Subsetting an object

Regular getindex (e.g. A[1:100, :]) and view work on all objects just as with an Array. view is always lazy, and reads from disk are deferred until getindex is used. DimensionalData.jl Dimensions and Selectors are the other way to subset an object, making use of the objects index to find values at e.g. certain X/Y coordinates. The available selectors are listed here:

At(x)get the index exactly matching the passed in value(s).
Near(x)get the closest index to the passed in value(s).
Where(f::Function)filter the array axis by a function of the dimension index values.
a..b/Between(a, b)get all indices between two values, excluding the high value.
Contains(x)get indices where the value x falls within an interval.

Use the .. selector to take a view of madagascar:

using Rasters, Plots
A = Raster(WorldClim{BioClim}, 5)
madagascar = view(A, X(43.25 .. 50.48), Y(-25.61 .. -12.04)) # Note the space between .. -12
plot(madagascar)

Methods that change the reslolution or extent of an object

Click through to the function documentation for more in-depth descriptions and examples.

aggregateaggregate data by the same or different amounts for each axis.
disaggregatesimilarly disaggregate data.
mosaicjoin rasters covering different extents into a single array or file.
cropshrink objects to specific dimension sizes or the extent of another object.
extendextend objects to specific dimension sizes or the extent of another object.
trimtrims areas of missing values for arrays and across stack layers.
resampleresample data to a different size and projection, or snap to another object.
warpuse gdalwarp on any object, e.g. a multidimensional NetCDF stack.

Methods that change an objects values:

Note that most regular Julia methods, such as replace, work as for a standard Array. These additional methods are commonly required in GIS applications.

classifyclassify values into categories.
maskmask an object by a polygon or Raster along X/Y, or other dimensions.
replace_missingreplace all missing values in an object and update missingval.

Point, polygon and table operation

rasterizerasterize points and geometries.
extractextract values from points or geometries.
inpolygonfind if a point or points are in a polygon.
zonalcalculate zonal statistics for an object masked by geometries.

Methods to load, write and modify data sources:

modifyreplace the data in objects. Useful to e.g. move objects to/from a GPU.
readread data to memory if it is on disk.
read!read data to predefined memory.
openopen the underlying data for manually reading or writing.
writewrite objects to file.

Altering and summarising arrays and stacks with regular julia methods

Most base methods work as for regular julia Arrays, such as reverse and rotations like rotl90. Base, statistics and linear algebra methods like mean that take a dims argument can also use the dimension name. To take the mean over the time dimension:

mean(A, dims=Ti)

broadcast works lazily from disk when lazy=true, and is only applied when data is directly indexed. Adding a dot to any function will use broadcast over a Raster just like an Array.

Broadcasting

For a disk-based array A, this will only be applied when indexing occurs or when we read the array.

A .*= 2

To broadcast directly to disk, we need to open the file in write mode:

open(Raster(filename); write=true) do O
    O .*= 2
end

To broadcast over a RasterStack use map, which applies a function to the raster layers of the stack.

newstack = map(stack) do raster
    raster .* 2
end

Modifying object properties

rebuild can be used to modify the fields of an object, generating a new object (but possibly holding the same arrays or files).

If you know that a file had an incorrectly specified missing value, you can do:

rebuild(A; missingval=-9999)

(replace_missing will actualy replace the current values)

Or if you need to change the name of the layer:

rebuild(A; name=:temperature)

set can be used to modify the nested properties of an objects dimensions, that are more difficult to change with rebuild. set works on the principal that dimension properties can only be in one specific field, so we generally don't have to specify which one it is. set will also try to update anything affected by a change you make.

This will set the X axis to specify points, instead of intervals:

set(A, X => Points)

We can also reassign dimensions, here X becomes Z:

set(A, X => Z)

setcrs(A, crs) and setmappedcrs(A, crs) will set the crs value/s of an object to any GeoFormat from GeoFormatTypes.jl.

Examples and Plotting

Plots.jl is fully supported by Rasters.jl, with recipes for plotting Raster and RasterStack provided. plot will plot a heatmap with axes matching dimension values. If mappedcrs is used, converted values will be shown on axes instead of the underlying crs values. contourf will similarly plot a filled contour plot.

Pixel resolution is limited to allow loading very large files quickly. max_res specifies the maximum pixel resolution to show on the longest axis of the array. It can be set manually to change the resolution (e.g. for large or high-quality plots):

using Rasters, Plots
A = Raster(WorldClim{BioClim}, 5)
plot(A; max_res=3000)

Loading and plotting data

Our first example simply loads a file from disk and plots it.

This netcdf file only has one layer, if it has more we could use RasterStack instead.

using Rasters, Plots
url = "https://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc";
filename = download(url, "tos_O1_2001-2002.nc");
A = Raster(filename)
180×170×24 Raster{Union{Missing, Float32},3} tos with dimensions: 
  X Mapped{Float64} Float64[1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Y Mapped{Float64} Float64[-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Ti Sampled{CFTime.DateTime360Day} CFTime.DateTime360Day[CFTime.DateTime360Day(2001-01-16T00:00:00), …, CFTime.DateTime360Day(2002-12-16T00:00:00)] ForwardOrdered Explicit Intervals
extent: Extent(X = (-0.0, 360.0), Y = (-80.0, 90.0), Ti = (CFTime.DateTime360Day(2001-01-01T00:00:00), CFTime.DateTime360Day(2003-01-01T00:00:00)))missingval: missing
crs: EPSG:4326
mappedcrs: EPSG:4326
values: [:, :, 1]
        -79.5       -78.5       …   86.5     87.5     88.5     89.5
   1.0     missing     missing     271.43   271.437  271.445  271.459
   3.0     missing     missing     271.431  271.438  271.445  271.459
   5.0     missing     missing     271.431  271.438  271.445  271.459
   7.0     missing     missing     271.431  271.439  271.446  271.459
   9.0     missing     missing  …  271.431  271.439  271.446  271.459
   ⋮                            ⋱                      ⋮      
 349.0     missing     missing     271.43   271.435  271.445  271.459
 351.0     missing     missing     271.43   271.435  271.445  271.459
 353.0     missing     missing     271.43   271.436  271.445  271.459
 355.0     missing     missing  …  271.43   271.436  271.445  271.459
 357.0     missing     missing     271.43   271.437  271.445  271.459
 359.0     missing     missing     271.43   271.437  271.445  271.459
[and 23 more slices...]

Objects with Dimensions other than X and Y will produce multi-pane plots. Here we plot every third month in the first year in one plot:

A[Ti=1:3:12] |> plot

Now plot the ocean temperatures around the Americas in the first month of 2001. Notice we are using lat/lon coordinates and date/time instead of regular indices. The time dimension uses DateTime360Day, so we need to load CFTime.jl to index it with Near.

using CFTime
A[Ti(Near(DateTime360Day(2001, 01, 17))), Y(-60.0 .. 90.0), X(45.0 .. 190.0)] |> plot

Now get the mean over the timespan, then save it to disk, and plot it as a filled contour:

Other plot functions and sliced objects that have only one X/Y/Z dimension fall back to generic DimensionalData.jl plotting, which will still correctly label plot axes.

using Statistics
# Take the mean
mean_tos = mean(A; dims=Ti)
180×170×1 Raster{Union{Missing, Float32},3} tos with dimensions: 
  X Mapped{Float64} Float64[1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Y Mapped{Float64} Float64[-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Ti Sampled{CFTime.DateTime360Day} CFTime.DateTime360Day(2002-01-16T00:00:00):Millisecond(2592000000):CFTime.DateTime360Day(2002-01-16T00:00:00) ForwardOrdered Explicit Intervals
extent: Extent(X = (-0.0, 360.0), Y = (-80.0, 90.0), Ti = (CFTime.DateTime360Day(2001-01-01T00:00:00), CFTime.DateTime360Day(2003-01-01T00:00:00)))missingval: missing
crs: EPSG:4326
mappedcrs: EPSG:4326
values: [:, :, 1]
        -79.5       -78.5       …   86.5     87.5     88.5     89.5
   1.0     missing     missing     271.427  271.434  271.443  271.454
   3.0     missing     missing     271.427  271.434  271.443  271.454
   5.0     missing     missing     271.427  271.434  271.443  271.454
   7.0     missing     missing     271.427  271.435  271.444  271.454
   9.0     missing     missing  …  271.427  271.435  271.444  271.454
   ⋮                            ⋱                      ⋮      
 349.0     missing     missing     271.427  271.433  271.443  271.454
 351.0     missing     missing     271.427  271.433  271.443  271.454
 353.0     missing     missing     271.427  271.433  271.443  271.454
 355.0     missing     missing  …  271.427  271.433  271.443  271.454
 357.0     missing     missing     271.427  271.433  271.443  271.454
 359.0     missing     missing     271.427  271.433  271.443  271.454

Plot a contour plot

contourf(mean_tos; dpi=300, size=(800, 400))

Write the mean values to disk

write("mean_tos.nc", mean_tos)
"mean_tos.nc"

Plotting recipes in DimensionalData.jl are the fallback for Rasters.jl when the object doesn't have 2 X/Y/Z dimensions, or a non-spatial plot command is used. So (as a random example) we could plot a transect of ocean surface temperature at 20 degree latitude :

A[Y(Near(20.0)), Ti(1)] |> plot

A basic species distribution modelling workflow

Load occurrences for the Mountain Pygmy Possum using GBIF.jl

using Rasters, GBIF2, Plots
records = GBIF2.occurrence_search("Burramys parvus"; limit=300)
300-element GBIF2.Table{GBIF2.Occurrence, JSON3.Array{JSON3.Object, Vector{UInt8}, SubArray{UInt64, 1, Vector{UInt64}, Tuple{UnitRange{Int64}}, true}}}┌──────────────────┬─────────────────┬─────────┬─────────┬─────────┬────────────
│ decimalLongitude │ decimalLatitude │    year │   month │     day │  kingdom  ⋯
│         Float64? │        Float64? │  Int64? │  Int64? │  Int64? │  String?  ⋯
├──────────────────┼─────────────────┼─────────┼─────────┼─────────┼────────────
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│          148.327 │        -36.2604 │    2016 │      11 │      15 │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│           148.22 │        -36.3511 │    2015 │      11 │      15 │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│          148.215 │        -36.4215 │    2011 │      11 │      21 │ Animalia  ⋯
│          missing │         missing │    2021 │       1 │       6 │ Animalia  ⋯
│          147.096 │        -36.9357 │    2020 │       2 │      10 │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│          148.347 │        -36.5047 │    2012 │      11 │      22 │ Animalia  ⋯
│          148.395 │        -36.5505 │    2016 │       1 │       3 │ Animalia  ⋯
│          148.236 │        -36.5249 │    2012 │      11 │      23 │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│          missing │         missing │ missing │ missing │ missing │ Animalia  ⋯
│        ⋮         │        ⋮        │    ⋮    │    ⋮    │    ⋮    │    ⋮      ⋱
└──────────────────┴─────────────────┴─────────┴─────────┴─────────┴────────────
                                                 77 columns and 285 rows omitted

Extract the longitude/latitude value to a Vector of points (a Tuple counts as a (x, y) point in GeoInterface.jl):

coords = [(r.decimalLongitude, r.decimalLatitude) for r in records if !ismissing(r.decimalLatitude)]
259-element Vector{Tuple{Float64, Float64}}:
 (148.326515, -36.260389)
 (148.2198, -36.351113)
 (148.214924, -36.421454)
 (147.096394, -36.935687)
 (148.347186, -36.504673)
 (148.395217, -36.550512)
 (148.235596, -36.524924)
 (148.240881, -36.400058)
 (151.250866, -33.831883)
 (147.1, -37.0)
 ⋮
 (148.378305, -36.413602)
 (148.391682, -36.373215)
 (148.368407, -36.413914)
 (148.307214, -36.452732)
 (148.256519, -36.442973)
 (148.461854, -36.009001)
 (148.461854, -36.009001)
 (148.461854, -36.009001)
 (148.40889, -36.006805)

Get BioClim layers and subset to south-east Australia

A = RasterStack(WorldClim{BioClim}, (1, 3, 7, 12))
se_aus = A[X(138 .. 155), Y(-40 .. -25), Band(1)]
RasterStack with dimensions: 
  X Projected{Float64} LinRange{Float64}(138.0, 154.833, 102) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(-25.1667, -39.8333, 89) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :bio1  Float32 dims: X, Y (102×89)
  :bio3  Float32 dims: X, Y (102×89)
  :bio7  Float32 dims: X, Y (102×89)
  :bio12 Float32 dims: X, Y (102×89)

Plot BioClim predictors and scatter occurrence points on all subplots

p = plot(se_aus);
kw = (legend=:none, opacity=0.5, markershape=:cross, markercolor=:black)
foreach(i -> scatter!(p, coords; subplot=i, kw...), 1:4)
display(p)

Then extract predictor variables and write to CSV.

using CSV
predictors = collect(extract(se_aus, coords))
CSV.write("burramys_parvus_predictors.csv", predictors)
"burramys_parvus_predictors.csv"

Or convert them to a DataFrame.

using DataFrames
df = DataFrame(predictors)
df[1:5, :]
5×5 DataFrame
Rowgeometrybio1bio3bio7bio12
Tuple…Float32Float32Float32Float32
1(148.327, -36.2604)8.8584641.039423.47881389.0
2(148.22, -36.3511)7.8357241.597523.50281500.0
3(148.215, -36.4215)7.8357241.597523.50281500.0
4(147.096, -36.9357)9.4083540.790523.08951292.0
5(148.347, -36.5047)8.420743.542523.1421223.0

Polygon masking, mosaic and plot

In this example we will mask the Scandinavian countries with border polygons, then mosaic together to make a single plot.

First, get the country boundary shape files using GADM.jl.

using Rasters, Shapefile, Plots, Dates, Downloads

# Download the shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "boundary_lines.shp"
Downloads.download(shapefile_url, shapefile_name)

# Load using Shapefile.jl
shapes = Shapefile.Handle(shapefile_name)
denmark_border = shapes.shapes[71]
norway_border = shapes.shapes[53]
sweden_border = shapes.shapes[54]
Polygon(4665 Points)

Then load raster data. We load some worldclim layers using RasterDataSources via Rasters.jl, and drop the Band dimension.

climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July)[Band(1)]
RasterStack with dimensions: 
  X Projected{Float64} LinRange{Float64}(-180.0, 179.833, 2160) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(89.8333, -90.0, 1080) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (2160×1080)
  :tmax Float32 dims: X, Y (2160×1080)
  :prec Int16 dims: X, Y (2160×1080)
  :wind Float32 dims: X, Y (2160×1080)

mask Denmark, Norway and Sweden from the global dataset using their border polygon, then trim the missing values. We pad trim with a 10 pixel margin.

mask_trim(climate, poly) = trim(mask(climate; with=poly); pad=10)

denmark = mask_trim(climate, denmark_border)
norway = mask_trim(climate, norway_border)
sweden = mask_trim(climate, sweden_border)
RasterStack with dimensions: 
  X Projected{Float64} LinRange{Float64}(9.5, 25.6667, 98) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(70.5, 53.6667, 102) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (98×102)
  :tmax Float32 dims: X, Y (98×102)
  :prec Int16 dims: X, Y (98×102)
  :wind Float32 dims: X, Y (98×102)

Plotting

First define a function to add borders to all subplots.

function borders!(p, poly)
    for i in 1:length(p)
        plot!(p, poly; subplot=i, fillalpha=0, linewidth=0.6)
    end
    return p
end
borders! (generic function with 1 method)

Now we can plot the individual countries.

dp = plot(denmark)
borders!(dp, denmark_border)
sp = plot(sweden)
borders!(sp, sweden_border)
np = plot(norway)
borders!(np, norway_border)

The Norway shape includes a lot of islands. Lets crop them out using .. intervals:

norway_region = climate[X(0..40), Y(55..73)]
plot(norway_region)

And mask it with the border again:

norway = mask_trim(norway_region, norway_border)
np = plot(norway)
borders!(np, norway_border)

Now we can combine the countries into a single raster using mosaic. first will take the first value if/when there is an overlap.

scandinavia = mosaic(first, denmark, norway, sweden)
RasterStack with dimensions: 
  X Projected{Float64} 3.1666666666666443:0.16666666666666666:32.49999999999998 ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} 72.66666666666666:-0.16666666666666666:52.99999999999999 ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (177×119)
  :tmax Float32 dims: X, Y (177×119)
  :prec Int16 dims: X, Y (177×119)
  :wind Float32 dims: X, Y (177×119)

And plot scandinavia, with all borders included:

p = plot(scandinavia)
borders!(p, denmark_border)
borders!(p, norway_border)
borders!(p, sweden_border)

And save to netcdf - a single multi-layered file, and tif, which will write a file for each stack layer.

write("scandinavia.nc", scandinavia)
write("scandinavia.tif", scandinavia)
("scandinavia_tmin.tif", "scandinavia_tmax.tif", "scandinavia_prec.tif", "scandinavia_wind.tif")

Rasters.jl provides a range of other methods that are being added to over time. Where applicable these methods read and write lazily to and from disk-based arrays of common raster file types. These methods also work for entire RasterStacks and RasterSeries using the same syntax.

Objects

Raster

Spatial raster data is essentially just an Array. But Raster wrappers allow treating them as an array that maintains its spatial index, crs and other metadata through all transformations. This means they can always be plotted and written to disk after applying most base Julia methods, and most broadcasts.

Rasters.AbstractRasterType
AbstractRaster <: DimensionalData.AbstractDimArray

Abstract supertype for objects that wrap an array (or location of an array) and metadata about its contents. It may be memory or hold a FileArray, which holds the filename, and is only opened when required.

AbstractRasters inherit from AbstractDimArray from DimensionalData.jl. They can be indexed as regular Julia arrays or with DimensionalData.jl Dimensions. They will plot as a heatmap in Plots.jl with correct coordinates and labels, even after slicing with getindex or view. getindex on a AbstractRaster will always return a memory-backed Raster.

source
Rasters.RasterType
Raster <: AbsractRaster

Raster(filepath::AbstractString, dims; kw...)
Raster(A::AbstractArray{T,N}, dims; kw...)
Raster(A::AbstractRaster; kw...)

A generic AbstractRaster for spatial/raster array data. It may hold memory-backed arrays or FileArray, that simply holds the String path to an unopened file. This will only be opened lazily when it is indexed with getindex or when read(A) is called. Broadcasting, taking a view, reversing and most other methods do not load data from disk: they are applied later, lazily.

Keywords

  • dims: Tuple of Dimensions for the array.
  • lazy: A Bool specifying if to load the stack lazily from disk. false by default.
  • name: Symbol name for the array, which will also retreive named layers if Raster is used on a multi-layered file like a NetCDF.
  • missingval: value reprsenting missing data, normally detected form the file. Set manually when you know the value is not specified or is incorrect. This will not change any values in the raster, it simply assigns which value is treated as missing. To replace all of the missing values in the raster, use replace_missing.
  • metadata: ArrayMetadata object for the array, or NoMetadata().
  • crs: the coordinate reference system of the objects XDim/YDim dimensions. Only set this if you know the detected crs is incrorrect, or it is not present in the file.
  • mappedcrs: the mapped coordinate reference system of the objects XDim/YDim dimensions. for Mapped lookups these are the actual values of the index. For Projected lookups this can be used to index in eg. EPSG(4326) lat/lon values, having it converted automatically. Only set this if the detected mappedcrs in incorrect, or the file does not have a mappedcrs, e.g. a tiff.

Internal Keywords

In some cases it is possible to set these keywords as well.

  • data: can replace the data in an AbstractRaster
  • refdims: Tuple of position Dimensions the array was sliced from, defaulting to ().
source
Rasters.RasterMethod
Raster(T::Type{<:RasterDataSource}, [layer]; kw...) => Raster

Load a RasterDataSource as a Raster. T and layers are passed to RasterDataSources.getraster, while kw args are for both getraster and Raster.

Keywords

  • month: an Int between 1 and 12, usually for Climate datasets
  • date: a DateTime object, usually for Weather datasets.
  • res: a String resolution, for datasets with multiple resolutions.

Other Raster keywords are passed to the Raster constructor.

See the docs for RasterDatasources.getraster for more specific details about data sources, layers and keyword arguments.

source

RasterStack

Spatial data often comes as a bundle of multiple named arrays, as in netcdf. RasterStack can represent this, or multiple files organised in a similar way.

Rasters.AbstractRasterStackType
AbstractRasterStack

Abstract supertype for objects that hold multiple AbstractRasters that share spatial dimensions.

They are NamedTuple-like structures that may either contain NamedTuple of AbstractRasters, string paths that will load AbstractRasters, or a single path that points to a file containing multiple layers, like NetCDF or HDF5. Use and syntax is similar or identical for all cases.

AbstractRasterStack can hold layers that share some or all of their dimensions. They cannot have the same dimension with different length or spatial extent as another layer.

getindex on an AbstractRasterStack generally returns a memory backed standard Raster. raster[:somelayer] |> plot plots the layers array, while raster[:somelayer, X(1:100), Band(2)] |> plot will plot the subset without loading the whole array.

getindex on an AbstractRasterStack with a key returns another stack with getindex applied to all the arrays in the stack.

source
Rasters.RasterStackType
RasterStack <: AbstrackRasterStack

RasterStack(data...; name, kw...)
RasterStack(data::Union{Vector,Tuple}; name, kw...)
RasterStack(data::NamedTuple; kw...))
RasterStack(s::AbstractRasterStack; kw...)
RasterStack(s::AbstractRaster; layersfrom=Band, kw...)
RasterStack(filename::AbstractString; kw...)

Load a file path or a NamedTuple of paths as a RasterStack, or convert arguments, a Vector or NamedTuple of Rasters to RasterStack.

Arguments

  • data: A NamedTuple of Rasters, or a Vector, Tuple or splatted arguments of Raster. The latter options must pass a name keyword argument.
  • filename: A file (such as netcdf or tif) to be loaded as a stack, or a directory path containing multiple files.

Keywords

  • name: Used as stack layer names when a Tuple, Vector or splat of Raster is passed in.
  • metadata: A Dict or DimensionalData.Metadata object.
  • refdims: Tuple of Dimension that the stack was sliced from.
  • layersfrom: Dimension to source stack layers from if the file is not already multi-layered. nothing is default, so that a single RasterStack(raster) is a single layered stack. RasterStack(raster; layersfrom=Band) will use the bands as layers.
  • lazy: A Bool specifying whether to load the stack lazily from disk. false by default.
files = (temp="temp.tif", pressure="pressure.tif", relhum="relhum.tif")
stack = RasterStack(files; mappedcrs=EPSG(4326))
stack[:relhum][Lat(Contains(-37), Lon(Contains(144))
source
Rasters.RasterStackMethod
RasterStack(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => RasterStack

Load a RasterDataSource as a RasterStack. T and layers are passed to RasterDataSources.getraster, while kw args are for both getraster and RasterStack.

Keywords

  • month: an Int between 1 and 12, usually for Climate datasets.
  • date: a DateTime object, usually for Weather datasets.
  • res: a String resolution, for datasets with multiple resolutions.

Other RasterStack keywords are passed to the RasterStack constructor.

See the docs for RasterDatasources.getraster for more specific details about data sources, layers and keyword arguments.

source

RasterSeries

A series is a meta-array that holds other files/data that is distributed over some dimension, often time. These files/data can be Rasters or RasterStacks.

Rasters.AbstractRasterSeriesType
AbstractRasterSeries <: DimensionalData.AbstractDimensionalArray

Abstract supertype for high-level DimensionalArray that hold RasterStacks, Rasters, or the paths they can be loaded from. RasterSeries are indexed with dimensions as with a AbstractRaster. This is useful when you have multiple files containing rasters or stacks of rasters spread over dimensions like time and elevation.

As much as possible, implementations should facilitate loading entire directories and detecting the dimensions from metadata.

This allows syntax like below for a series of stacks of arrays:

RasterSeries[Time(Near(DateTime(2001, 1))][:temp][Y(Between(70, 150)), X(Between(-20,20))] |> plot`

RasterSeries is the concrete implementation.

source
Rasters.RasterSeriesType
RasterSeries <: AbstractRasterSeries

RasterSeries(arrays::AbstractArray{<:AbstractRaster}, dims; kw...)
RasterSeries(stacks::AbstractArray{<:AbstractRasterStack}, dims; kw...)
RasterSeries(filepaths::AbstractArray{<:AbstractString}, dims; child, duplicate_first, kw...)
RasterSeries(dirpath:::AbstractString, dims; ext, child, duplicate_first, kw...)

Concrete implementation of AbstractRasterSeries. A RasterSeries is an array of Rasters or RasterStacks, along some dimension(s).

Arguments

  • dims: series dimension/s.

Keywords

  • refdims: existing reference dimension/s.
  • child: constructor of child objects for use when filenames are passed in, can be Raster or RasterStack. Defaults to Raster.
  • lazy: load files lazily. This is true by default for series, as it is common to load many files to openare over lazily.
  • duplicate_first::Bool: wether to duplicate the dimensions and metadata of the first file with all other files. This can save load time with a large series where dimensions are essentially identical. false by default.
  • ext: filename extension such as ".tiff" or ".nc". Use if only a directory path is passed in.
  • kw: keywords passed to the child constructor Raster or RasterStack if only file names are passed in.
source
Rasters.RasterSeriesMethod
RasterSeries(T::Type{<:RasterDataSource}, [layers::Union{Symbol,AbstractArray,Tuple}]; kw...) => AbstractRasterSeries

Load a RasterDataSource as an AbstractRasterSeries. T, args are passed to RasterDataSource.getraster, while kw args are for both getraster and RasterSeries.

Keywords

  • month: a Vector or range of Int between 1 and 12, usually for Climate datasets.
  • date: a Vector of DateTime objects, usually for Weather datasets.
  • res: a String resolution, for datasets with multiple resolutions.

Other RasterSeries keywords are passed to the RasterSeries constructor.

See the docs for RasterDatasources.getraster for more specific details about data sources, layers and keyword arguments.

source

Dimensions

Rasters uses X, Y, and Z dimensions from DimensionalData.jl to represent spatial directions like longitude, latitude and the vertical dimension, and subset data with them. Ti is used for time, and Band represent bands. Other dimensions can have arbitrary names, but will be treated generically. See DimensionalData.jl for more details on how they work.

Rasters.BandType
Band <: Dimension

Band(val=:)

Band Dimension for multi-band rasters.

Example:

banddim = Band(10:10:100)
# Or
val = A[Band(1)]
# Or
mean(A; dims=Band)
source

Lookup Arrays

These specify properties of the index associated with e.g. the X and Y dimension. Rasters.jl defines additional lookup arrays: Projected to handle dimensions with projections, and Mapped where the projection is mapped to another projection like EPSG(4326). Mapped is largely designed to handle NetCDF dimensions, especially with Explicit spans.

Rasters.ProjectedType
Projected <: AbstractProjected

Projected(order, span, sampling, crs, mappedcrs)
Projected(; order=AutoOrder(), span=AutoSpan(), sampling=AutoSampling(), crs, mappedcrs=nothing)

An AbstractSampled LookupArray with projections attached.

Fields and behaviours are identical to Sampled with the addition of crs and mappedcrs fields.

If both crs and mappedcrs fields contain CRS data (in a GeoFormat wrapper from GeoFormatTypes.jl) the selector inputs and plot axes will be converted from and to the specified mappedcrs projection automatically. A common use case would be to pass mappedcrs=EPSG(4326) to the constructor when loading eg. a GDALarray:

GDALarray(filename; mappedcrs=EPSG(4326))

The underlying crs will be detected by GDAL.

If mappedcrs is not supplied (ie. mappedcrs=nothing), the base index will be shown on plots, and selectors will need to use whatever format it is in.

source
Rasters.MappedType
Mapped <: AbstractProjected

Mapped(order, span, sampling, crs, mappedcrs)
Mapped(; order=AutoOrder(), span=AutoSpan(), sampling=AutoSampling(), crs=nothing, mappedcrs)

An AbstractSampled LookupArray, where the dimension index has been mapped to another projection, usually lat/lon or EPSG(4326). Mapped matches the dimension format commonly used in netcdf files.

Fields and behaviours are identical to Sampled with the addition of crs and mappedcrs fields.

The mapped dimension index will be used as for Sampled, but to save in another format the underlying crs may be used to convert it.

source

Data sources

Rasters.jl uses a number of backends to load raster data. Raster, RasterStack and RasterSeries will detect which backend to use for you, automatically.

GRD

R GRD files can be loaded natively, using Julias MMap - which means they are very fast, but are not compressed. They are always 3 dimensional, and have Y, X and Band dimensions.

NetCDF

NetCDF .nc files are loaded using NCDatasets.jl. Layers from files can be loaded as Raster("filename.nc"; key=:layername). Without key the first layer is used. RasterStack("filename.nc") will use all netcdf variables in the file that are not dimensions as layers.

NetCDF layers can have arbitrary dimensions. Known, common dimension names are converted to X, Y Z, and Ti, otherwise Dim{:layername} is used. Layers in the same file may also have different dimensions.

NetCDF files still have issues loading directly from disk for some operations. Using read(ncstack) may help.

GDAL

All files GDAL can access, such as .tiff and .asc files, can be loaded, using ArchGDAL.jl. These are generally best loaded as Raster("filename.tif"), but can be loaded as RasterStack("filename.tif"; layersfrom=Band), taking layers from the Band dimension, which is also the default.

SMAP

The Soil Moisture Active-Passive dataset provides global layers of soil moisture, temperature and other related data, in a custom HDF5 format. Layers are always 2 dimensional, with Y and X dimensions.

These can be loaded as multi-layered RasterStack("filename.h5"). Individual layers can be loaded as Raster("filename.h5"; key=:layerkey), without key the first layer is used.

Rasters.smapseriesFunction
smapseries(filenames::AbstractString; kw...)
smapseries(filenames::Vector{<:AbstractString}, dims=nothing; kw...)

RasterSeries loader for SMAP files and whole folders of files, organised along the time dimension. Returns a RasterSeries.

Arguments

  • filenames: A String path to a directory of SMAP files, or a vector of String paths to specific files.
  • dims: Tuple containing Ti dimension for the series. Automatically generated form filenames unless passed in.

Keywords

  • kw: Passed to RasterSeries.
source

Writing file formats to disk

Files can be written to disk in all formats other than SMAP HDF5 using write("filename.ext", A). See the docs for write. They can (with some caveats) be written to different formats than they were loaded in as, providing file-type conversion for spatial data.

Some metadata may be lost in formats that store little metadata, or where metadata conversion has not been completely implemented.

RasterDataSources.jl integration

RasterDataSources.jl standardises the download of common raster data sources, with a focus on datasets used in ecology and the environmental sciences. RasterDataSources.jl is tightly integrated into Rasters.jl, so that datsets and keywords can be used directly to download and load data as a Raster, RasterStack, or RasterSeries.

using Rasters, Plots, Dates
A = Raster(WorldClim{Climate}, :tavg; month=June)
plot(A)

See the docs for Raster, RasterStack and RasterSeries, and the docs for RasterDataSources.getraster for syntax to specify various data sources.

Exported functions

Rasters.jl is a direct extension of DimensionalData.jl. See DimensionalData.jl docs for the majority of types and functions that can be used in Rasters.jl.

Functions more specific to geospatial data are included in Rasters.jl, and listed below.

Rasters.aggregateFunction
aggregate(method, object, scale; filename, progress, skipmissing)

Aggregate a Raster, or all arrays in a RasterStack or RasterSeries, by scale using method.

Arguments

  • method: a function such as mean or sum that can combine the value of multiple cells to generate the aggregated cell, or a Locus like Start() or Center() that specifies where to sample from in the interval.
  • object: Object to aggregate, like AbstractRasterSeries, AbstractStack, AbstractRaster or Dimension.
  • scale: the aggregation factor, which can be an integer, a tuple of integers for each dimension, or any Dimension, Selector or Int combination you can usually use in getindex. Using a Selector will determine the scale by the distance from the start of the index.

When the aggregation scale of is larger than the array axis, the length of the axis is used.

Keywords

  • progress: show a progress bar.
  • skipmissingval: if true, any missingval will be skipped during aggregation, so that only areas of all missing values will be aggregated to missingval. If false, any aggegrated area containing a missingval will be assigned missingval.
  • filename: a filename to write to directly, useful for large files.
  • suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(st) are the default.

Example

using Rasters, Statistics, Plots
using Rasters: Center
st = read(RasterStack(WorldClim{Climate}; month=1))
ag = aggregate(Center(), st, (Y(20), X(20)); skipmissingval=true, progress=false)
plot(ag)
savefig("build/aggregate_example.png"); nothing
# output

aggregate

Note: currently it is faster to aggregate over memory-backed arrays. Use read on src before use where required.

source
Rasters.aggregate!Function
aggregate!(method, dst::AbstractRaster, src::AbstractRaster, scale; skipmissingval=false)

Aggregate array src to array dst by scale, using method.

Arguments

  • method: a function such as mean or sum that can combine the value of multiple cells to generate the aggregated cell, or a Locus like Start() or Center() that species where to sample from in the interval.
  • scale: the aggregation factor, which can be an integer, a tuple of integers for each dimension, or any Dimension, Selector or Int combination you can usually use in getindex. Using a Selector will determine the scale by the distance from the start of the index in the src array.

When the aggregation scale of is larger than the array axis, the length of the axis is used.

Keywords

  • progress: show a progress bar.
  • skipmissingval: if true, any missingval will be skipped during aggregation, so that only areas of all missing values will be aggregated to missingval. If false, any aggegrated area containing a missingval will be assigned missingval.

Note: currently it is much faster to aggregate over memory-backed arrays. Use read on src before use where required.

source
Rasters.boolmaskFunction
boolmask(obj::Raster; [missingval])
boolmask(obj; [to, res, size])

Create a mask array of Bool values, from another Raster. An AbstractRasterStack or AbstractRasterSeries are also accepted, but a mask is taken of the first layer or object not all of them.

The array returned from calling boolmask on a AbstractRaster is a Raster with the same dimensions as the original array and a missingval of false.

Arguments

  • obj: a Raster, a GeoInterface.jl geometry, or a vector or table of geometries.

Raster / RasterStack Keywords

  • missingval: The missing value of the source array, with default missingval(raster).

Geometry keywords

  • to: a Raster, RasterStack, Tuple of Dimension or Extents.Extent. If no to object is provided the extent will be calculated from the geometries, Additionally, when no to object or an Extent is passed for to, the size or res keyword must also be used.

  • res: the resolution of the dimensions, a Real or Tuple{<:Real,<:Real}.

  • size: the size of the output array, as a Tuple{Int,Int} or single Int for a square. Only required when to is not used or is anExtents.Extent, otherwisesize`.

  • shape: Force data to be treated as :polygon, :line or :point geometries. using points or lines as polygons may have unexpected results.

  • boundary: for polygons, include pixels where the :center is inside the polygon, where the line :touches the pixel, or that are completely :inside inside the polygon. The default is :center.

And specifically for shape=:polygon:

  • boundary: include pixels where the :center is inside the polygon, where the line :touches the pixel, or that are completely :inside inside the polygon. The default is :center.

Example

using Rasters, Plots, Dates
wc = Raster(WorldClim{Climate}, :prec; month=1)
boolmask(wc) |> plot

savefig("build/boolmask_example.png"); nothing

# output

boolmask

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.classifyFunction
classify(x, pairs; lower=(>=), upper=(<), others=nothing)
classify(x, pairs...; lower, upper, others)

Create a new array with values in x classified by the values in pairs.

pairs can hold tuples fo values (2, 3), a Fix2 function e.g. <=(1), a Tuple of Fix2 e.g. (>=(4), <(7)), or an IntervalSets.jl interval, e.g. 3..9 or OpenInterval(10, 12). pairs can also be a n * 3 matrix where each row is lower bounds, upper bounds, replacement.

If tuples or a Matrix are used, the lower and upper keywords define how the lower and upper boundaries are chosen.

If others is set other values not covered in pairs will be set to that values.

Arguments

  • x: a Raster or RasterStack
  • pairs: each pair contains a value and a replacement, a tuple of lower and upper range and a replacement, or a Tuple of Fix2 like (>(x), <(y).

Keywords

  • lower: Which comparison (< or <=) to use for lower values, if Fix2 are not used.
  • upper: Which comparison (> or >=) to use for upper values, if Fix2 are not used.
  • others: A value to assign to all values not included in pairs. Passing nothing (the default) will leave them unchanged.

Example

using Rasters, Plots
A = Raster(WorldClim{Climate}, :tavg; month=1)
classes = <=(15) => 10,
          15..25 => 20,
          25..35 => 30,
          >(35) => 40
classified = classify(A, classes; others=0, missingval=0)
plot(classified; c=:magma)

savefig("build/classify_example.png"); nothing

# output

classify

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.classify!Function
classify!(x, pairs...; lower, upper, others)
classify!(x, pairs; lower, upper, others)

Classify the values of x in-place, by the values in pairs.

If Fix2 is not used, the lower and upper keywords

If others is set other values not covered in pairs will be set to that values.

Arguments

  • x: a Raster or RasterStack
  • pairs: each pair contains a value and a replacement, a tuple of lower and upper range and a replacement, or a Tuple of Fix2 like (>(x), <(y).

Keywords

  • lower: Which comparison (< or <=) to use for lower values, if Fix2 are not used.
  • upper: Which comparison (> or >=) to use for upper values, if Fix2 are not used.
  • others: A value to assign to all values not included in pairs. Passing nothing (the default) will leave them unchanged.

Example

classify! to disk, with key steps:

  • copying a tempory file so we don't write over the RasterDataSources.jl version.
  • use open with write=true to open the file with disk-write permissions.
  • use Float32 like 10.0f0 for all our replacement values and other, because the file is stored as Float32. Attempting to write some other type will fail.
using Rasters, Plots, RasterDataSources
# Download and copy the file
filename = getraster(WorldClim{Climate}, :tavg; month=6)
tempfile = tempname() * ".tif"
cp(filename, tempfile)
# Define classes
classes = (5, 15) => 10,
          (15, 25) => 20,
          (25, 35) => 30,
          >=(35) => 40
# Open the file with write permission
open(Raster(tempfile); write=true) do A
    classify!(A, classes; others=0)
end
# Open it again to plot the changes
plot(Raster(tempfile); c=:magma)

savefig("build/classify_bang_example.png"); nothing

# output

classify!

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.convertlookupFunction
convertlookup(dstlookup::Type{<:LookupArray}, x)

Convert the dimension lookup between Projected and Mapped. Other dimension lookups pass through unchanged.

This is used to e.g. save a netcdf file to GeoTiff.

source
Rasters.cropFunction
crop(x; to)
crop(xs...; to)

Crop one or multiple AbstractRaster or AbstractRasterStack x to match the size of the object to, or smallest of any dimensions that are shared.

Otherwise crop to the size of the keyword argument to. This can be a Tuple of Dimension or any object that will return one from dims(to).

crop is lazy, using a view into the object rather than alocating new memory.

Keywords

  • to: the object to crop to. If to keyword is passed, the smallest shared area of all xs is used.

As crop is lazy, filename and suffix keywords don't apply.

Example

Crop to another raster:

using Rasters, Plots
evenness = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
rnge = Raster(EarthEnv{HabitatHeterogeneity}, :range)

# Roughly cut out New Zealand from the evenness raster
nz_bounds = X(165 .. 180), Y(-50 .. -32)
nz_evenness = evenness[nz_bounds...]

# Crop range to match evenness
nz_range = crop(rnge; to=nz_evenness)
plot(nz_range)

savefig("build/nz_crop_example.png")
nothing

# output

new zealand evenness cropped

Crop to a polygon:

using Rasters, Plots, Dates, Shapefile, Downloads

# Download a borders shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "boundary.shp"
isfile(shapefile_name) || Downloads.download(shapefile_url, shapefile_name)
shp = Shapefile.Handle(shapefile_name).shapes[6]

evenness = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
argentina_evenness = crop(evenness; to=shp)
plot(argentina_evenness)

savefig("build/argentina_crop_example.png"); nothing

# output

argentina evenness cropped

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.crsFunction
crs(x)

Get the projected coordinate reference system of a Y or X Dimension, or of the Y/X dims of an AbstractRaster.

For Mapped lookup this may be nothing as there may be no projected coordinate reference system at all.

source
Rasters.disaggregateFunction
disaggregate(method, object, scale; filename, progress, keys)

Disaggregate array, or all arrays in a stack or series, by some scale.

Arguments

  • method: a function such as mean or sum that can combine the value of multiple cells to generate the aggregated cell, or a Locus like Start() or Center() that species where to sample from in the interval.
  • object: Object to aggregate, like AbstractRasterSeries, AbstractStack, AbstractRaster or a Dimension.
  • scale: the aggregation factor, which can be an integer, a tuple of integers for each dimension, or any Dimension, Selector or Int combination you can usually use in getindex. Using a Selector will determine the scale by the distance from the start of the index.

Keywords

  • progress: show a progress bar.

Note: currently it is faster to aggregate over memory-backed arrays. Use read on src before use where required.

source
Rasters.disaggregate!Function
disaggregate!(method, dst::AbstractRaster, src::AbstractRaster, filename, scale)

Disaggregate array src to array dst by some scale, using method.

  • method: a function such as mean or sum that can combine the value of multiple cells to generate the aggregated cell, or a Locus like Start() or Center() that species where to sample from in the interval.
  • scale: the aggregation factor, which can be an integer, a tuple of integers for each dimension, or any Dimension, Selector or Int combination you can usually use in getindex. Using a Selector will determine the scale by the distance from the start of the index in the src array.

Note: currently it is faster to aggregate over memory-backed arrays. Use read on src before use where required.

source
Rasters.extendFunction
extend(xs...; [to])
extend(xs; [to])
extend(x::Union{AbstractRaster,AbstractRasterStack}; to, kw...)

Extend one or multiple AbstractRaster to match the area covered by all xs, or by the keyword argument to.

Keywords

  • to: the Raster or dims to extend to. If no to keyword is passed, the largest shared area of all xs is used.
  • atol: the absolute tolerance value to use when comparing the index of x and to.
  • filename: a filename to write to directly, useful for large files.
  • suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(st) are the default.
using Rasters, Plots
evenness = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
rnge = Raster(EarthEnv{HabitatHeterogeneity}, :range)

# Roughly cut out South America
sa_bounds = X(-88 .. -32), Y(-57 .. 13)
sa_evenness = evenness[sa_bounds...]

# Extend range to match the whole-world raster
sa_range = extend(sa_evenness; to=rnge)
plot(sa_range)

savefig("build/extend_example.png")
nothing
# output

extend

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.extractFunction
extract(x, geoms; atol)

Extracts the value of Raster or RasterStack at given points, returning an iterable of NamedTuple with properties for :geometry and raster or stack layer values.

Note that if objects have more dimensions than the length of the point tuples, sliced arrays or stacks will be returned instead of single values.

Arguments

  • x: a Raster or RasterStack to extract values from.
  • geoms: GeoInterface.jl compatible geometries, or tables or iterables of geometries.

Keywords

  • atol: a tolorerance for floating point lookup values for when the LookupArray contains Points. atol is ignored for Intervals.

Example

Here we extact points matching the occurrence of the Mountain Pygmy Possum, Burramis parvus. This could be used to fit a species distribution model.

using Rasters, GBIF2, CSV

# Get a stack of BioClim layers, and replace missing values with `missing`
st = RasterStack(WorldClim{BioClim}, (1, 3, 5, 7, 12))[Band(1)] |> replace_missing

# Download some occurrence data
obs = GBIF2.occurrence_search("Burramys parvus"; limit = 10)

# Convert observations to points
pnts = collect((o.decimalLongitude, o.decimalLatitude) for o in obs if !ismissing(o.decimalLongitude))

# use `extract` to get values for all layers at each observation point.
# We `collect` to get a `Vector` from the lazy iterator.
collect(extract(st, pnts))

# output
5-element Vector{NamedTuple{(:geometry, :bio1, :bio3, :bio5, :bio7, :bio12), Tuple{Tuple{Float64, Float64}, Float32, Float32, Float32, Float32, Float32}}}:
 (geometry = (148.326515, -36.260389), bio1 = 8.8584585, bio3 = 41.039413, bio5 = 21.6955, bio7 = 23.47875, bio12 = 1389.0)
 (geometry = (148.2198, -36.351113), bio1 = 7.8357186, bio3 = 41.597527, bio5 = 20.51325, bio7 = 23.50275, bio12 = 1500.0)
 (geometry = (148.214924, -36.421454), bio1 = 7.8357186, bio3 = 41.597527, bio5 = 20.51325, bio7 = 23.50275, bio12 = 1500.0)
 (geometry = (147.096394, -36.935687), bio1 = 9.408354, bio3 = 40.790546, bio5 = 22.39425, bio7 = 23.0895, bio12 = 1292.0)
 (geometry = (148.347186, -36.504673), bio1 = 8.420698, bio3 = 43.542496, bio5 = 21.20175, bio7 = 23.142, bio12 = 1223.0)
source
Rasters.inpolygonFunction
inpolygon(points, poly)

Check if a point or points are inside a polygon.

This algorithm is efficient for many points, less so for a single point.

Arguments

  • points: a GeoInterface.jl compatable object, or AbstractVector of objects.
  • poly: a GeoInterface.jl compatable polygon, or a feature or feature collection of polygons.

Returns a Bool or AbstractVector{Bool}.

source
Rasters.mappedcrsFunction
mappedcrs(x)

Get the mapped coordinate reference system for the Y/X dims of an array.

In Projected lookup this is used to convert Selector values form the mappedcrs defined projection to the underlying projection, and to show plot axes in the mapped projection.

In Mapped lookup this is the coordinate reference system of the index values.

source
Rasters.mappedindexFunction
mappedindex(x)

Get the index value of a dimension converted to the mappedcrs value.

Whithout ArchGDAL loaded, this is just the regular dim value.

source
Rasters.maskFunction
mask(A:AbstractRaster; with, missingval=missingval(A))
mask(x; with)

Return a new array with values of A masked by the missing values of with, or by the shape of with, if with is a geometric object.

Arguments

  • x: a Raster or RasterStack

Keywords

  • with: an AbstractRaster, or any GeoInterface.jl compatible objects or table. The coordinate reference system of the point must match crs(A).
  • missingval: the missing value to use in the returned file.
  • filename: a filename to write to directly, useful for large files.
  • suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(st) are the default.

Geometry keywords

These can be used when with is a GeoInterface.jl compatible object:

  • shape: Force data to be treated as :polygon, :line or :point geometries. using points or lines as polygons may have unexpected results.
  • boundary: for polygons, include pixels where the :center is inside the polygon, where the line :touches the pixel, or that are completely :inside inside the polygon. The default is :center.

Example

Mask an unmasked AWAP layer with a masked WorldClim layer, by first resampling the mask.

using Rasters, Plots, Dates

# Load and plot the file
awap = read(Raster(AWAP, :tmax; date=DateTime(2001, 1, 1)))
a = plot(awap; clims=(10, 45))

# Create a mask my resampling a worldclim file
wc = Raster(WorldClim{Climate}, :prec; month=1)
wc_mask = resample(wc; to=awap)

# Mask
awap_masked = mask(awap; with=wc_mask)
b = plot(awap_masked; clims=(10, 45))

savefig(a, "build/mask_example_before.png");
savefig(b, "build/mask_example_after.png"); nothing
# output

Before mask:

before mask

After mask:

after mask

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.mask!Function
mask!(x; with, missingval=missingval(A))

Mask A by the missing values of with, or by all values outside with if it is a polygon.

If with is a polygon, creates a new array where points falling outside the polygon have been replaced by missingval(A).

Return a new array with values of A masked by the missing values of with, or by a polygon.

Arguments

  • x: a Raster or RasterStack.

Keywords

  • with: another AbstractRaster, a AbstractVector of Tuple points, or any GeoInterface.jl AbstractGeometry. The coordinate reference system of the point must match crs(A).
  • missingval: the missing value to write to A in masked areas, by default missingval(A).

Example

Mask an unmasked AWAP layer with a masked WorldClim layer, by first resampling the mask to match the size and projection.

using Rasters, Plots, Dates

# Load and plot the file
awap = read(RasterStack(AWAP, (:tmin, :tmax); date=DateTime(2001, 1, 1)))
a = plot(awap; clims=(10, 45), c=:imola)

# Create a mask my resampling a worldclim file
wc = Raster(WorldClim{Climate}, :prec; month=1)
wc_mask = resample(wc; to=awap)

# Mask
mask!(awap; with=wc_mask)
b = plot(awap; clims=(10, 45))

savefig(a, "build/mask_bang_example_before.png");
savefig(b, "build/mask_bang_example_after.png"); nothing

# output

Before mask!:

before mask!

After mask!:

after mask!

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.missingmaskFunction
missingmask(obj::Raster; kw...)
missingmask(obj; [to, res, size])

Create a mask array of missing and true values, from another Raster. AbstractRasterStack or AbstractRasterSeries are also accepted, but a mask is taken of the first layer or object not all of them.

For AbstractRaster the default missingval is missingval(A), but others can be chosen manually.

The array returned from calling missingmask on a AbstractRaster is a Raster with the same size and fields as the original array.

Keywords

  • to: a Raster, RasterStack, Tuple of Dimension or Extents.Extent. If no to object is provided the extent will be calculated from the geometries, Additionally, when no to object or an Extent is passed for to, the size or res keyword must also be used.

  • res: the resolution of the dimensions, a Real or Tuple{<:Real,<:Real}.

  • size: the size of the output array, as a Tuple{Int,Int} or single Int for a square. Only required when to is not used or is anExtents.Extent, otherwisesize`.

  • shape: Force data to be treated as :polygon, :line or :point geometries. using points or lines as polygons may have unexpected results.

  • boundary: for polygons, include pixels where the :center is inside the polygon, where the line :touches the pixel, or that are completely :inside inside the polygon. The default is :center.

Example

using Rasters, Plots, Dates
wc = Raster(WorldClim{Climate}, :prec; month=1)
missingmask(wc) |> plot

savefig("build/missingmask_example.png"); nothing

# output

missingmask

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.mosaicFunction
mosaic(f, regions...; missingval, atol)
mosaic(f, regions; missingval, atol)

Combine regions into a single raster.

Arguments

  • f: A reducing function (mean, sum, first, last etc.) for values where regions overlap.
  • regions: Iterable or splatted Raster or RasterStack.

Keywords

  • missingval: Fills empty areas, and defualts to the missingval of the first region.
  • atol: Absolute tolerance for comparison between index values. This is often required due to minor differences in range values due to floating point error. It is not applied to non-float dimensions. A tuple of tolerances may be passed, matching the dimension order.
  • filename: a file to write to directly.

If your mosaic has has apparent line errors, increase the atol value.

Example

Here we cut out Australia and Africa from a stack, and join them with mosaic.

using Rasters, Plots
st = RasterStack(WorldClim{Climate}; month=1);

africa = st[X(-20.0 .. 60.0), Y(-40.0 .. 35.0)]
a = plot(africa)

aus = st[X(100.0 .. 160.0), Y(-50.0 .. -10.0)]
b = plot(aus)

# Combine with mosaic
mos = mosaic(first, aus, africa)
c = plot(mos)

savefig(a, "build/mosaic_example_africa.png")
savefig(b, "build/mosaic_example_aus.png")
savefig(c, "build/mosaic_example_combined.png")
nothing
# output

Individual continents

arica

aus

Mosaic of continents

mosaic

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.mosaic!Function
mosaic!(f, x, regions...; missingval, atol)
mosaic!(f, x, regions::Tuple; missingval, atol)

Combine regions in x using the function f.

Arguments

  • f a function (e.g. mean, sum, first or last) that is applied to values where regions overlap.
  • x: A Raster or RasterStack. May be a an opened disk-based Raster, the result will be written to disk. With the current algorithm, the read speed is slow.
  • regions: source objects to be joined. These should be memory-backed (use read first), or may experience poor performance. If all objects have the same extent, mosaic is simply a merge.

Keywords

  • missingval: Fills empty areas, and defualts to the `missingval/ of the first layer.
  • atol: Absolute tolerance for comparison between index values. This is often required due to minor differences in range values due to floating point error. It is not applied to non-float dimensions. A tuple of tolerances may be passed, matching the dimension order.

Example

Cut out Australia and Africa stacks, then combined them into a single stack.

using Rasters, Statistics, Plots
st = read(RasterStack(WorldClim{Climate}; month=1))
aus = st[X=100.0 .. 160.0, Y=-50.0 .. -10.0]
africa = st[X=-20.0 .. 60.0, Y=-40.0 .. 35.0]
mosaic!(first, st, aus, africa)
plot(st)
savefig("build/mosaic_bang_example.png")
nothing
# output

mosaic

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.pointsFunction
points(A::AbstractRaster; dims=(YDim, XDim), ignore_missing) => Array{Tuple}

Returns a generator of the points in A for dimensions in dims, where points are a tuple of the values in each specified dimension index.

Keywords

  • dims the dimensions to return points from. The first slice of other layers will be used.
  • ignore_missing: wether to ignore missing values in the array when considering points. If true, all points in the dimensions will be returned, if false only the points that are not === missingval(A) will be returned.

The order of dims determines the order of the points.

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.rasterizeFunction
rasterize(obj; to, fill, kw...)

Rasterize the a GeoInterface.jl compatable geometry or feature, or a Tables.jl table with a :geometry column of GeoInterface.jl objects, or X, Y points columns.

Arguments

  • obj: a GeoInterface.jl AbstractGeometry, or a nested Vector of AbstractGeometry, or a Tables.jl compatible object containing points and values columns.

Keywords

These are detected automatically from obj where possible.

  • fill: the value to fill a polygon with. A Symbol or tuple of Symbol will be used to retrieve properties from features or column values from table rows.
  • atol: an absolute tolerance for rasterizing to dimensions with Points sampling.
  • filename: a filename to write to directly, useful for large files.
  • suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(st) are the default.

Geometry keywords

These can be used when a GeoInterface.AbstractGeometry is passed in.

  • to: a Raster, RasterStack, Tuple of Dimension or Extents.Extent. If no to object is provided the extent will be calculated from the geometries, Additionally, when no to object or an Extent is passed for to, the size or res keyword must also be used.

  • res: the resolution of the dimensions, a Real or Tuple{<:Real,<:Real}.

  • size: the size of the output array, as a Tuple{Int,Int} or single Int for a square. Only required when to is not used or is anExtents.Extent, otherwisesize`.

  • shape: Force data to be treated as :polygon, :line or :point geometries. using points or lines as polygons may have unexpected results.

  • boundary: for polygons, include pixels where the :center is inside the polygon, where the line :touches the pixel, or that are completely :inside inside the polygon. The default is :center.

Example

Rasterize a shapefile for China and plot, with a border.

using Rasters, Plots, Dates, Shapefile, Downloads
using Rasters.LookupArrays

# Download a borders shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "country_borders.shp"
isfile(shapefile_name) || Downloads.download(shapefile_url, shapefile_name)

# Loade the shapes for china
china_border = Shapefile.Handle(shapefile_name).shapes[10]

# Rasterize the border polygon
china = rasterize(china_border; res=0.1, missingval=0, fill=1, boundary=:touches)

# And plot
p = plot(china; color=:spring)
plot!(p, china_border; fillalpha=0, linewidth=0.6)

savefig("build/china_rasterized.png"); nothing

# output

rasterize

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.rasterize!Function
rasterize!(dest, data; fill, atol)

Rasterize the geometries in data into the Raster or RasterStack x, using the values specified by fill.

Arguments

  • dest: a Raster or RasterStack to rasterize into.
  • data: an GeoInterface.jl compatible object or AbstractVector of such objects, or a Tables.jl compatible table containing GeoInterface compatible objects or columns with point names X and Y.
  • fill: the value to fill a polygon with. A Symbol or tuple of Symbol will be used to retrieve properties from features or column values from table rows.

Keywords

These are detected automatically from A and data where possible.

  • atol: an absolute tolerance for rasterizing points to dimensions with Points sampling.

Geometry keywords

These can be used when a GeoInterface.AbstractGeometry is passed in.

  • shape: Force data to be treated as :polygon, :line or :point.

And specifically for shape=:polygon:

  • boundary: include pixels where the :center is inside the polygon, where the line :touches the pixel, or that are completely :inside inside the polygon.

Table keywords

  • name: A Symbol to return a Raster from a single column, or Tuple of Symbol to return a RasterStack from multiple columns.

Example

using Rasters, Plots, Dates, Shapefile, GeoInterface, Downloads
using Rasters.LookupArrays

# Download a borders shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "country_borders.shp"
isfile(shapefile_name) || Downloads.download(shapefile_url, shapefile_name)

# Load the shapes for indonesia
indonesia_border = Shapefile.Handle(shapefile_name).shapes[1]

# Make an empty EPSG 4326 projected Raster of the area of Indonesia
dimz = X(90.0:0.1:145; mode=Projected(; sampling=Points(), crs=EPSG(4326))),
       Y(-15.0:0.1:10.9; mode=Projected(; sampling=Points(), crs=EPSG(4326)))

A = Raster(zeros(UInt32, dimz); missingval=UInt32(0))

# Rasterize each island with a different number
for (i, shp) in enumerate(GeoInterface.getring(indonesia_border))
    rasterize!(A, shp; fill=i)
end

# And plot
p = plot(Rasters.trim(A); color=:spring)
plot!(p, indonesia_border; fillalpha=0, linewidth=0.7)

savefig("build/indonesia_rasterized.png"); nothing

# output

rasterize

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.resampleFunction
resample(x, resolution::Number; crs, method)
resample(x; to, method)
resample(xs...; to=first(xs), method)

resample uses ArchGDAL.gdalwarp to resample a Raster or RasterStack to a new resolution and optionally new crs, or to snap to the bounds, resolution and crs of the object to.

Arguments

  • x: the object to resample.
  • resolution: a Number specifying the resolution for the output. If the keyword argument crs (described below) is specified, resolution must be in units of the crs.

Keywords

  • to: an AbstractRaster whose resolution, crs and bounds will be snapped to. For best results it should roughly cover the same extent, or a subset of A.

  • crs: A GeoFormatTypes.GeoFormat such as EPSG(x) or WellKnownText(string) specifying an output crs (A will be reprojected to crs in addition to being resampled). Defaults to crs(A)

  • method: A Symbol or String specifying the method to use for resampling. From the docs for gdalwarp:

    • :near: nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).
    • :bilinear: bilinear resampling.
    • :cubic: cubic resampling.
    • :cubicspline: cubic spline resampling.
    • :lanczos: Lanczos windowed sinc resampling.
    • :average: average resampling, computes the weighted average of all non-NODATA contributing pixels. rms root mean square / quadratic mean of all non-NODATA contributing pixels (GDAL >= 3.3)
    • :mode: mode resampling, selects the value which appears most often of all the sampled points.
    • :max: maximum resampling, selects the maximum value from all non-NODATA contributing pixels.
    • :min: minimum resampling, selects the minimum value from all non-NODATA contributing pixels.
    • :med: median resampling, selects the median value of all non-NODATA contributing pixels.
    • :q1: first quartile resampling, selects the first quartile value of all non-NODATA contributing pixels.
    • :q3: third quartile resampling, selects the third quartile value of all non-NODATA contributing pixels.
    • :sum: compute the weighted sum of all non-NODATA contributing pixels (since GDAL 3.1)

    Where NODATA values are set to missingval.

Note: missingval of missing does not work with GDAL. Use replace_missing(A, newmissingval) to assign a missing value before using resample if the current value is missing. This will be automated in future versions.

Example

Resample a WorldClim layer to match an EarthEnv layer:

using Rasters, Plots
A = Raster(WorldClim{Climate}, :prec; month=1)
B = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)

a = plot(A)
b = plot(resample(A; to=B))

savefig(a, "build/resample_example_before.png");
savefig(b, "build/resample_example_after.png"); nothing
# output

Before resample:

before resample

After resample:

after resample

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.replace_missingFunction
replace_missing(a::AbstractRaster, newmissingval)
replace_missing(a::AbstractRasterStack, newmissingval)

Replace missing values in the array or stack with a new missing value, also updating the missingval field/s.

Keywords

  • filename: a filename to write to directly, useful for large files.
  • suffix: a string or value to append to the filename. A tuple of suffix will be applied to stack layers. keys(st) are the default.

Example

using Rasters
A = Raster(WorldClim{Climate}, :prec; month=1) |> replace_missing
missingval(A)
# output
missing
source
Rasters.reprojectFunction
reproject(target::GeoFormat, x)

Reproject the dimensions of x to a different crs.

Arguments

  • target: any crs in a GeoFormatTypes.jl wrapper, e.g. EPSG, WellKnownText, ProjString.
  • x: a Dimension, Tuple of Dimension, Raster or RasterStack.

Dimensions without an AbstractProjected lookup (such as a Ti dimension) are silently returned without modification.

source
reproject(source::GeoFormat, target::GeoFormat, dim::Dimension, val)

reproject uses ArchGDAL.reproject, but implemented for a reprojecting a value array of values, a single dimension at a time.

source
Rasters.setcrsFunction
setcrs(x, crs)

Set the crs of a Raster, RasterStack, Tuple of Dimension, or a Dimension.

source
Rasters.setmappedcrsFunction
setmappedcrs(x, crs)

Set the mapped crs of a Raster, a RasterStack, a Tuple of Dimension, or a Dimension.

source
Base.skipmissingFunction
skipmissing(itr::Raster)

Returns an iterable over the elements in a Raster object, skipping any values equal to either the missingval or missing.

source
Rasters.sliceFunction
slice(A::Union{AbstractRaster,AbstractRasterStack,AbstracRasterSeries}, dims) => RasterSeries

Slice an object along some dimension/s, lazily using view.

For a single Raster or RasterStack this will return a RasterSeries of Raster or RasterStack that are slices along the specified dimensions.

For a RasterSeries, the output is another series where the child objects are sliced and the series dimensions index is now of the child dimensions combined. slice on a RasterSeries with no dimensions will slice along the dimensions shared by both the series and child object.

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.subsetFunction
subset(s::AbstractRasterStack, keys)

Subset a stack to hold only the layers in keys, where keys is a Tuple or Array of String or Symbol, or a Tuple or Array of Int

Depreciated: As it is now possible to subset NamedTuple by indexing with a Tuple of Symbol in getindex, that is also possible for any AbstractDimStack like RasterStack. So subset is obsolete and will be remove in future versions.

Use:

```julia s[(:key1, :key2)]

source
Rasters.trimFunction
trim(x; dims::Tuple, pad::Int)

Trim missingval(x) from x for axes in dims, returning a view of x.

Arguments

  • x: A Raster or RasterStack. For stacks, all layers must having missing values for a pixel for it to be trimmed.

Keywords

  • dims: By default dims=(XDim, YDim), so that trimming keeps the area of X and Y that contains non-missing values along all other dimensions.
  • pad: The trimmed size will be padded by pad on all sides, although padding will not be added beyond the original extent of the array.

trim does not accept filename/suffix arguments as it does not alter the underlying data.

Example

Create trimmed layers of Australian habitat heterogeneity.

using Rasters, Plots
layers = (:evenness, :range, :contrast, :correlation)
st = RasterStack(EarthEnv{HabitatHeterogeneity}, layers)

# Roughly cut out australia
ausbounds = X(100 .. 160), Y(-50 .. -10)
aus = st[ausbounds...]
a = plot(aus)

# Trim missing values and plot
b = plot(trim(aus))

savefig(a, "build/trim_example_before.png");
savefig(b, "build/trim_example_after.png"); nothing

# output

Before trim:

before trim

After trim:

after trim

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.warpFunction
warp(A::AbstractRaster, flags::Dict)

Gives access to the GDALs gdalwarp method given a Dict of flag => value arguments that can be converted to strings, or vectors where multiple space-separated arguments are required.

Arrays with additional dimensions not handled by GDAL (other than X, Y, Band) are sliced, warped, and then combined to match the original array dimensions.

See the gdalwarp docs for a list of arguments.

Example

This simply resamples the array with the :tr (output file resolution) and :r flags, giving us a pixelated version:

using Rasters, RasterDataSources, Plots
A = Raster(WorldClim{Climate}, :prec; month=1)
a = plot(A)

flags = Dict(
    :tr => [2.0, 2.0],
    :r => :near,
)
b = plot(warp(A, flags))

savefig(a, "build/warp_example_before.png");
savefig(b, "build/warp_example_after.png"); nothing

# output

Before warp:

before warp

After warp:

after warp

In practise, prefer resample for this. But warp may be more flexible.

WARNING: This feature is experimental. It may change in future versions, and may not be 100% reliable in all cases. Please file github issues if problems occur.

source
Rasters.zonalFunction
zonal(f, x::Union{Raster,RasterStack}; of, kw...)

Calculate zonal statistics for the the zone of a Raster or RasterStack covered by the of object/s.

Arguments

  • f: any function that reduces an iterable to a single value, such as sum or Statistics.mean
  • x: A Raster or RasterStack
  • of: A Raster, RasterStack, dim tuple, extent, GeoInterface.jl compatible geometry, Tables.jl compatible table of a :geometry column, or an AbstractVector of any of these objects..

Keywords

These can be used when of is a GeoInterface.jl compatible object:

  • shape: Force data to be treated as :polygon, :line or :point, where possible.
  • boundary: for polygons, include pixels where the :center is inside the polygon, where the line :touches the pixel, or that are completely :inside inside the polygon. The default is :center.

Example

using Rasters, Shapefile, DataFrames, Downloads, Statistics, Dates

# Download a borders shapefile
ne_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries"
shp_url, dbf_url  = ne_url * ".shp", ne_url * ".dbf"
shp_name, dbf_name = "country_borders.shp", "country_borders.dbf"
isfile(shp_name) || Downloads.download(shp_url, shp_name)
isfile(dbf_url) || Downloads.download(dbf_url, dbf_name)

# Download and read a raster stack from WorldClim
st = RasterStack(WorldClim{Climate}; month=Jan, lazy=false)

# Load the shapes for world countries
countries = Shapefile.Table(shp_name) |> DataFrame

# Calculate the january mean of all climate variables for all countries
january_stats = zonal(mean, st; of=countries, boundary=:touches) |> DataFrame

# Add the country name column (natural earth has some string errors it seems)
insertcols!(january_stats, 1, :country => first.(split.(countries.ADMIN, r"[^A-Za-z ]")))

# output
258×8 DataFrame
 Row │ country                       tmin           tmax          tavg         ⋯
     │ SubStrin…                     Float32?       Float32?      Float32?     ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │ Indonesia                          21.539         29.1784       25.3587 ⋯
   2 │ Malaysia                           21.2955        28.4185       24.8569
   3 │ Chile                               7.24828       17.9385       12.5933
   4 │ Bolivia                            17.233         27.7662       22.4995
   5 │ Peru                               15.0116        25.5386       20.2751 ⋯
   6 │ Argentina                          13.683         27.6834       20.6831
   7 │ Dhekelia Sovereign Base Area  missing        missing       missing
   8 │ Cyprus                              5.03622       14.2194        9.6276
  ⋮  │              ⋮                      ⋮             ⋮              ⋮      ⋱
 252 │ Spratly Islands                    25.1           29.4          27.2    ⋯
 253 │ Clipperton Island             missing        missing       missing
 254 │ Macao S                       missing        missing       missing
 255 │ Ashmore and Cartier Islands   missing        missing       missing
 256 │ Bajo Nuevo Bank               missing        missing       missing      ⋯
 257 │ Serranilla Bank               missing        missing       missing
 258 │ Scarborough Reef              missing        missing       missing
                                                  5 columns and 243 rows omitted
source

File operations

These Base and DimensionalData methods have specific Rasters.jl versions:

DimensionalData.modifyFunction
modify(f, series::AbstractRasterSeries)

Apply function f to the data of the child object. If the child is an AbstractRasterStack the function will be passed on to its child AbstractRasters.

f must return an idenically sized array.

This method triggers a complete rebuild of all objects, and disk based objects will be transferred to memory.

An example of the usefulnesss of this is for swapping out array backend for an entire series to CuArray from CUDA.jl to copy data to a GPU.

source
Base.openFunction
open(f, A::AbstractRaster; write=false)

open is used to open any lazy=true AbstractRaster and do multiple operations on it in a safe way. The write keyword opens the file in write lookup so that it can be altered on disk using e.g. a broadcast.

f is a method that accepts a single argument - an Raster object which is just an AbstractRaster that holds an open disk-based object. Often it will be a do block:

lazy=false (in-memory) rasters will ignore open and pass themselves to f.

# A is an `Raster` wrapping the opened disk-based object.
open(Raster(filepath); write=true) do A
    mask!(A; with=maskfile)
    A[I...] .*= 2
    # ...  other things you need to do with the open file
end

By using a do block to open files we ensure they are always closed again after we finish working with them.

source
Base.readFunction
read(A::AbstractRaster)
read(A::AbstractRasterStack)
read(A::AbstractRasterSeries)

read will move a Rasters.jl object completely to memory.

source
Base.read!Function
read!(src::Union{AbstractString,AbstractRaster}, dst::AbstractRaster)
read!(src::Union{AbstractString,AbstractRasterStack}, dst::AbstractRasterStack)
read!(scr::AbstractRasterSeries, dst::AbstractRasterSeries)

read! will copy the data from src to the object dst.

src can be an object or a file-path String.

source
Base.writeFunction
Base.write(filename::AbstractString, A::AbstractRaster; kw...)

Write an AbstractRaster to file, guessing the backend from the file extension.

Keyword arguments are passed to the write method for the backend.

source
Base.write(filename::AbstractString, s::AbstractRasterStack; suffix, kw...)

Write any AbstractRasterStack to file, guessing the backend from the file extension.

Keywords

  • suffix: suffix to append to file names. By default the layer key is used.

Other keyword arguments are passed to the write method for the backend.

If the source can't be saved as a stack-like object, individual array layers will be saved.

source
Base.write(filename::AbstractString, ::Type{GRDfile}, s::AbstractRaster)

Write a Raster to a .grd file with a .gri header file. The extension of filename will be ignored.

Returns filename.

source
Base.write(filename::AbstractString, ::Type{NCDfile}, A::AbstractRaster)

Write an NCDarray to a NetCDF file using NCDatasets.jl

Returns filename.

source
Base.write(filename::AbstractString, ::Type{NCDfile}, s::AbstractRasterStack; kw...)

Write an NCDstack to a single netcdf file, using NCDatasets.jl.

Currently Metadata is not handled for dimensions, and Metadata from other AbstractRaster @types is ignored.

Keywords

Keywords are passed to NCDatasets.defVar.

  • append: If true, the variable of the current Raster will be appended to filename. Note that the variable of the current Raster should be not exist before. If not, you need to set append = false. Rasters.jl can not overwrite a previous existing variable.
  • fillvalue: A value filled in the NetCDF file to indicate missing data. It will be stored in the _FillValue attribute.
  • chunksizes: Vector integers setting the chunk size. The total size of a chunk must be less than 4 GiB.
  • deflatelevel: Compression level: 0 (default) means no compression and 9 means maximum compression. Each chunk will be compressed individually.
  • shuffle: If true, the shuffle filter is activated which can improve the compression ratio.
  • checksum: The checksum method can be :fletcher32 or :nochecksum (checksumming is disabled, which is the default)
  • typename (string): The name of the NetCDF type required for vlen arrays (https://web.archive.org/save/https://www.unidata.ucar.edu/software/netcdf/netcdf-4/newdocs/netcdf-c/nc005fdef005fvlen.html)
source
Base.write(filename::AbstractString, ::Type{GDALfile}, A::AbstractRaster; kw...)

Write a Raster to file using GDAL.

Keywords

  • driver: A GDAL driver name or a GDAL driver retrieved via ArchGDAL.getdriver(drivername). Guessed from the filename extension by default.

  • options::Dict{String,String}: A dictionary containing the dataset creation options passed to the driver. For example: Dict("COMPRESS"=>"DEFLATE")

    Valid options for the drivers can be looked up here: https://gdal.org/drivers/raster/index.html

Returns filename.

source

Internals

Rasters.FileArrayType
FileArray{X} <: DiskArrays.AbstractDiskArray

Filearray is a DiskArrays.jl AbstractDiskArray. Instead of holding an open object, it just holds a filename string that is opened lazily when it needs to be read.

source
Rasters.FileStackType
FileStack{X,K}

FileStack{X,K}(filename, types, sizes, eachchunk, haschunks, write)

A wrapper object that holds file pointer and size/chunking metadata for a multi-layered stack stored in a single file, typically netcdf or hdf5.

X is a backend type like NCDfile, and K is a tuple of Symbol keys.

source
Rasters.OpenStackType
OpenStack{X,K}

OpenStack{X,K}(dataset)

A wrapper for any stack-like opened dataset that can be indexed with Symbol keys to retrieve AbstractArray layers.

OpenStack is usually hidden from users, wrapped in a regular RasterStack passed as the function argument in open(stack) when the stack is contained in a single file.

X is a backend type like NCDfile, and K is a tuple of Symbol keys.

source
Rasters.RasterDiskArrayType
RasterDiskArray <: DiskArrays.AbstractDiskArray

A basic DiskArrays.jl wrapper for objects that don't have one defined yet. When we open a FileArray it is replaced with a RasterDiskArray.

source